{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "[![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=notebooks/89_image_array_viz.ipynb)\n",
    "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/notebooks/89_image_array_viz.ipynb)\n",
    "[![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)\n",
    "\n",
    "**Visualizing in-memory raster datasets and image arrays**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install \"leafmap[raster]\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import leafmap\n",
    "import rasterio\n",
    "import rioxarray\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "Download two sample raster datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "url1 = \"https://opengeos.org/data/raster/landsat.tif\"\n",
    "url2 = \"https://opengeos.org/data/raster/srtm90.tif\"\n",
    "satellite = leafmap.download_file(url1, \"landsat.tif\", overwrite=True)\n",
    "dem = leafmap.download_file(url2, \"srtm90.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "The Landsat image contains 3 bands: nir, red, and green. Let's calculate NDVI using the nir and red bands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = rasterio.open(satellite)\n",
    "nir = dataset.read(4).astype(float)\n",
    "red = dataset.read(1).astype(float)\n",
    "ndvi = (nir - red) / (nir + red)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Create an in-memory raster dataset from the NDVI array and use the projection and extent of the Landsat image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_image = leafmap.array_to_image(ndvi, source=satellite)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "Visualize the Landsat image and the NDVI image on the same map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_raster(satellite, indexes=[4, 1, 2], vmin=0, vmax=120, layer_name=\"Landsat 7\")\n",
    "m.add_raster(ndvi_image, colormap=\"Greens\", layer_name=\"NDVI\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "You can also specify the image metadata (e.g., cellsize, crs, and transform) when creating the in-memory raster dataset.\n",
    "\n",
    "First, check the metadata of the origina image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.profile"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "Check the crs of the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Check the transform of the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.transform"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "Create an in-memory raster dataset from the NDVI array and specify the cellsize, crs, and transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = (30.0, 0.0, -13651650.0, 0.0, -30.0, 4576290.0)\n",
    "ndvi_image = leafmap.array_to_image(\n",
    "    ndvi, cellsize=30, crs=\"EPSG:3857\", transform=transform\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "Add the NDVI image to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_raster(satellite, indexes=[4, 1, 2], vmin=0, vmax=120, layer_name=\"Landsat 7\")\n",
    "m.add_raster(ndvi_image, colormap=\"Greens\", layer_name=\"NDVI\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "Use rioxarray to read raster datasets into xarray DataArrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = rioxarray.open_rasterio(dem)\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "Classify the DEM into 2 elevation classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "array = ds.sel(band=1)\n",
    "masked_array = xr.where(array < 2000, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25",
   "metadata": {},
   "source": [
    "Visualize the DEM and the elevation class image on the same map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_raster(dem, colormap=\"terrain\", layer_name=\"DEM\")\n",
    "m.add_raster(masked_array, colormap=\"coolwarm\", layer_name=\"Classified DEM\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27",
   "metadata": {},
   "source": [
    "Add a split map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[37.6, -119], zoom=9)\n",
    "m.split_map(\n",
    "    dem,\n",
    "    masked_array,\n",
    "    left_args={\n",
    "        \"layer_name\": \"DEM\",\n",
    "        \"colormap\": \"terrain\",\n",
    "    },\n",
    "    right_args={\n",
    "        \"layer_name\": \"Classified DEM\",\n",
    "        \"colormap\": \"coolwarm\",\n",
    "    },\n",
    ")\n",
    "m"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
